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ABSTRACT 



The optimal minimum lime control (i.e. bang-bang controller) is applied to 
the fast reaction missile defense problem. From Pontryagin, the optimal control 
was determined to be a function of the adjoint in the minimization of the 
Hamiltonian. The control may also be posed either as a function of time or as a 
function of the states. The state space can be partitioned into regions, surfaces 
and curves where the optimal control action is either its maximum plus or minus 
N. 

In missile simulation problems, the method of adjoints is often used in 
parametric studies of errors and miss distance. This technique is developed 
both graphically and mathematically, and is used here to help one visualize the 
solution trajectory and families of optimal trajectories for all possible initial 
conditions. 
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I. INTRODUCTION 



In this era of technological advancement, more and more effort is being 
devoted to automation to increase speed and efficiency. Everything is 
becoming faster, smaller, more efficient and more effective, from fuzzy logic 
controlled appliances to smart weapon systems. As our engineering and design 
of physical systems improves we are able to generate faster and more accurate 
mechanical devices, and the control systems for such devices must improve as 
well. In pushing to develop the fastest, most efficient controller for whatever 
our application, we introduce the bang-bang controller. 

Many of our current defensive missile systems were designed to be used 
against threats that are no longer of greatest importance. A point defense 
missile may now be required to bring down a target that has a significant speed 
advantage, and be able to do it in such a way as to control the fallout after the 
impact. 

Even if we redesign our systems, current economic conditions make it 
unlikely that we could build up an arsenal of new weapon systems. However, 
we could upgrade our currently existing arsenal by replacing the original control 
logic with a newer, more capable controller. 

In this report we will develop second and third order minimum time optimal 
controllers and apply them to a fast reaction missile defense problem. While we 
are using the missile defense problem as our example, it should be noted that 
the controllers are not limited to this example and may be applied as generic 
controllers for a wide variety of systems. 
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II. SECOND ORDER CONTROLLER 



A. MINIMUM TIME OPTIMAL CONTROL 

Designing a control system is frequently a hit-and-miss process using a 
variety of design techniques to iteratively create a system that meets specific 
criteria. The performance of the system may be defined in the time and 
frequency domain in terms of overshoot, settling time, etc., or it may be defined 
by some externally measured criteria. For example, a satellite control system 
may be designed to minimize fuel consumption. 

"The objective of optimal control theory is to determine the control signals 
that will cause a process to satisfy the physical constraints and at the same 
time minimize (or maximize) some performance criterion." [1] 

1. Problem Definition 

The system and optimization problem is defined as 

x = Ax-f-Bu (2.1) 

with x(0) = c, and x(tf) = 0, while minimizing 



"I 

J = Jdt. 



( 2 . 2 ) 



Let us consider the simple example with 



X = 



0 1 
0 0 



x-f- 



(2.3) 



From Pontryagin [1], we find we can minimize J by minimizing the Hamiltonian 

H = l-l-pjX2 -t-p2U. (2.4) 
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This is minimum when u is operating at its maximum possible value and with 
opposite the sign of the adjoint p 2 . Thus we have 

u = -N-sign(p 2 ) (2.5) 



where 



This has a solution 



P 



dx 





Pi =Ci 

P2 =-Clt + C2. 

A typical solution would appear as seen in Fig. 2.1. 



( 2 . 6 ) 



(2.7) 

( 2 . 8 ) 



P2’U 




Figure 2.1 Solution to Minimized Hamiltonian 

Note that from the adjoint solution, the control may change sign only 
once. We would like to solve this problem for all possible initial conditions and 
only one terminal condition (x(tf) = 0). Hence it makes sense to look at this 
problem in negative time (adjoint), starting from the end point of motion. From 
the uniqueness theorem in ordinary differential equations, only two trajectories 
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can emanate from the origin; one for u = -N and one for u = +N, as shown in 
Figure 2.2. In this second order example, our solution is constrained to the Xi, 
X 2 plane. These negative time trajectories divide the state space (here a plane) 
into two parts. Solution trajectories emanating from these curves constitute all 
possible solutions for all possible initial conditions. 




Figure 2.2 Zero Trajectory Curves 

2. Minimum Time Trajectories, Second Order 
Consider the system 





"0 


r 




"O' 


X = 


0 


0 


x + 


1 



This can be represented in flow diagram form as 



(2.9) 



u 

* 



1 I 

S Xj s 

> • > • 
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The control, u, may be defined for a bang-bang control system as ±N, where N 
is some constant. Given any starting values for the states, there are only two 
possible paths for the states to take; one corresponding to u = -i-N, and one 
corresponding to u = -N. If we desire to drive all states to zero, and we know 
that we can only apply u = ±N, Figure 2.3 shows the paths followed by the 
states given various starting values. 





X2 

1 1 'W 


(v 


V / i 



Figure 2.3 Minimum Time Trajectory Solution Curves 



3. Second Order Solution Trajectories 




This system has a solution 




x(t) = ^(t,0)x(0) + A(t,0)u(0) 


(2.10) 


where 






(2.11) 



and 
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( 2 . 12 ) 



A= 

Jo 



"1 O' 




'0 r 


1 




'0 1 




'0 r 


0 




+ 




t+— 










r+... 


_0 1_ 




O 

O 


2! 


O 

O 




O 

O 




'1 O' 




'0 t' 


1 




0 0' 






'1 t' 




+ 




H 






t 


^+...= 




0 1 




0 0 


2! 




0 0 






0 1 



Expanding e'^‘, we get 



e^‘ = 



1 + At + — y 
2! 



(2.13) 



Since is the zero matrix, all higher orders of A will also be the zero matrix 
and we may drop them. Evaluating A (for u fixed) we find 



A = Bdi = [' 

Jo Jo 



'1 t-x' 


'0' 


_0 1 


1 



dx 



't-x' 




'tX-^T^' 


t 




1 


dx = 




= 


2 




X 


0 


t 



The solution for this system is 



(2.14) 



x(t) = 



1 t 
0 1 



x(0) + 






u(0), 



(2.15) 



or in scalar equations 

x,(t) = Xi(0) + X2(0)t + ^u(0)t^ (2.16) 

X2(t) = X2(0) + u(0)t. (2.17) 

These equations describe the states as a function of time given any initial 
conditions and the fixed control effort, u. 



B. ANALYTIC SOLUTION FOR SECOND ORDER SWITCHING 
TIME 

We may now treat this system as a boundary value problem and 
analytically solve for switching times of the control effort. 
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1. Solving Boundary Value Problems 

Since the control is piecewise constant, (±N), we can separate the 
problem into two pieces and match boundary values at the point where u 
changes sign. In other words, if the system moves from point x(to) through 
point x(ts) to point x(tf), we can solve the problem in two parts. Each of our 
boundary value problems can be stated in such a way as to supply simplifying 
boundary conditions, i.e., setting initial or final values to zero. Optimal control 
theory tells us that for a second order system there will be at most one 
switching time, the change from u = +N to u = -N, or visa versa. Let us 
consider first the solution for our system with some initial condition [X](0), 
X2(0)], as shown in Figure 2.4. 




Figure 2.4 Minimum Time Trajectory From a Fixed Initial 
Condition 

The time from t=0 to ts is the length of time before the control effort 
changes sign. We cannot solve directly for ts because we do not have enough 
information on the boundary values. Since we are solving this problem for 
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arbitrary initial conditions, xi(0) and X2(0), xi(ts) and X 2 (ts) are unknown. We 
do know, however, that the final states, xi(tf) and X 2 (tf), will be zero, and from 
this we can determine the zero-trajectory curve in negative time from the 
origin. This curve, shown in Figure 2.5., is the only curve that will go through 
the origin with u = +N. Therefore the switching time will occur when the path 
of a u = -N curve intersects this curve. There are an infinite number of u = -N 
curves intersecting this curve, however, only one curve will go through any 
particular set of initial conditions. 

To simplify the problem, we will start with initial conditions on the Xi 
axis. Define Xi(0) = ci, X2(0) = 0, xi(tf) = X 2 (tf) = 0, and u = -N. This situation 
is described in Figure 2.6 where ti may be defined as the time from to to tj. 
Equation (2.16) may now be written 

Xj(t) = Xi(0)-l-X2(0)t-l- jut^ = Xj(0)-^Nt^ (2.18) 



Since the curves for u = -i-N and u = -N are symmetric, it can be noted that 





Xi(ti) = ^X|(0) 


(2.19) 


and therefore 








ix,(0)=x,(0)-iNt,2 


(2.20) 


and finally 







ti = 




( 2 . 21 ) 



This is valid for any initial conditions such that xi(0) > 0, X2(0) = 0 and u = -N. 
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Figure 2.5 Simplified Boundary Value Problem #1 




Figure 2.6 Simplified Boundary Value Problem #2 
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The next simplification involves getting from the positive X 2 axis down 
to the positive Xi axis as shown in Figure 2.7. Given Xi(0) = 0, X2(0) = positive 
real, xi(tf) = positive real, X 2 (tf) = 0, and u = -N. 

Equation (2.17) may be written 

X2(t) = X2(0)-Nt. (2.22) 



At t = t 2 , X 2 (t 2 ) = 0 giving 

0 = X2(0)-Nt2 



(2.23) 



and 



_X2(0) 

N • 



(2.24) 



The final stage is to translate the initial condition off of the X 2 axis to 
some unknown initial condition. The time required for the states to get from 
some initial condition X 2 ( 0 ) to X 2 (t) is not based on the xi state, and is therefore 
always equal to t 2 , as seen in Figure 2.8. 

We solve (2.16) for Xi(t 2 ) 



Xl(t2) =Xi(0)+X2(0)t2-yNt2 



= Xi(0) + X2(0)<'''2^°^'' 



N 7 



-In 



X2(0) J 



= Xi(0) + — X2(0). 
IV 2N ^ 



(2.25) 
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Figure 2.7 Simplified Boundary Value Problem #3 




Figure 2.8 Simplified Boundary Value Problem #4 
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Finally we combine these solutions to determine the switching time, tj, 
for any initial conditions such that Xi(0) > 0 and X2(0) > 0 with u = -N. From 



Figure 2.3 we see that xi(0) for the t] solution is the same point as Xi(t 2 ) from 
the t 2 solution so that 



t. = 



and 





'x.(t,) 

N 


x,(0)+ x^(0) 

i 2N 


\ N 


k(0) 1 

V N 2 


f X2(0)V 
1 N j 



(2.26) 



ts = t,+t2 = 



|x,(0) ^ lfx,(0) 



x,(0) 

N 



(2.27) 



, N 2\ N J 
These equations were derived for specific initial conditions and are 
only valid where the initial conditions lie above the zero trajectory curves of 
Figure 2.2. In order to expand the capabilities of the control system for all initial 
conditions we re-derive the switching time equations for initial conditions 
Xj(0) < 0, X2(0) < 0, and u(0) = +N for the solution 



l-x,(0) 1 




' x,(0) 


y N 2 


1 N j 


N 



(2.28) 



This equation is valid only when the initial values are below the zero trajectory 
curves of Figure 2.2. For this discussion we will limit ourselves to initial 
conditions above the zero trajectory curves and use the solution (2.27). 
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2. Simulation of Analytic Switching Time Controller 

Using a computer simulation, we test the accuracy of the solution by 
choosing initial conditions and observing the response to our control effort. We 
simulated the example system of (2.9) with xi(0) > 0, X2(0) > 0, and N = -1. 
The control effort, u, changes sign to -N at the calculated switching time, tj. 
The results show that the states pass through the origin of the state space and 
the error at the origin is associated with the discretization of the simulation, as 
shown in Figures 2.9 and 2.10. The switching time shown in Figure 2.10 is the 
calculated value. The control effort actually switches at the first sampling time 
following the caluculated switching time, ts. If we increase the sample rate for 
the simulation we reduce the terminal miss error. Upon reaching the origin 
some additional control logic must be devised to maintain position. 

While the switching time solution is a minimum time solution for driving 
the states to the origin, or any desired values, it must be shut off at tf. The 
switching time solution cannot adapt to a time varying situation as the switching 
time is defined solely on the initial conditions. 
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XI 



Figure 2.9 Simulation of Analytic Switching Time Solution 
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Figure 2.10 Control Effort for Analytic Switching Time 
Simulation 
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C. SOLUTION TO SECOND ORDER SWITCHING CURVES 



An option for improving the capabilities of the second order controller is to 
remove the time dependancy of the control effort. The control problem may be 
posed as a function of the states instead of a function of time for some simple 
systems. By defining the control effort as only a function of the states the 
control effort is updated continuously and may change immediately in response 
to a change in the system parameters. 

There are only two possible control efforts, -N and +N. Therefore, any 
position in space will have either +N or -N control effort to drive it along its 
unique path to the origin, and we can solve this system for the second order 
switching curves that divide the state space into two parts, see Figure 2.2. 

1. Second Order Switching Laws 
Given the system 



’X]' 




’0 


r 


’^r 


+ 


"O' 


>2. 




_0 


o_ 


.^2. 




_1_ 



we can integrate the state equations. Setting u = ±N 



Choosing the boundary values so that Xi(tf) = X 2 (tf) = 0 we may rewrite the 
equations as 




Xi(tf) = ±|Ntf +X2(0)tf + x,( 0 ) 
X2(tf) = ±Ntf+X2(0) 



(2.32) 

(2.33) 



or 



0 — iyNtf + X2(0)tf + Xj(0) 



(2.34) 
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Solving for tf 



0 = ±Ntf + X2(0). 



_^X2(0) 



tf = + 



N 



Then for any value of t 



t = + ^ . 

N 



Substituting this into equation (2.34) 

0 = Xi(t) + X2(t) 







':pX2(t)^ 


1 N J 


2 1 


. N j 



^ N ^ N 



(2.35) 

(2.36) 

(2.37) 

(2.38) 

(2.39) 



2 

0 = x,(t) + |^^. (2.40) 

We make the substitution 

+ x|(t) = -X 2 (t)|x 2 (t)| (2.41) 

so that 

0 = Xi(t)-^X2(t)|x2(t)|. (2.42) 

Equation (2.42) describes the two parabolic trajectories referred to as 
the zero trajectory curves because any states on these curves will, with the 
appropriately signed control effort, be driven to the origin. These zero 
trajectory curves divide the state space into two regions of opposite control 
effort. If the states are above the curves then u = -N, but if they are below the 
curves then u = +N. Therefore our control effort, u, may be defined 

u =-Nsign(xi(t)-^X 2 (t)|x 2 (t)|). (2.43) 



16 



With this switching law we drive the states from any initial conditions to the 
origin with at most one change of the control effort. 

For example, let us define the initial conditions of the states to be 
above the zero trajectory curves so that Xi(0) < 0 and X 2 ( 0 ) > 0, as shown in 
Figure 2.11. From (2.43) the control effort is u = -N, which drives the states 
along the parabolic trajectory shown in Figure 2.11. 

When the states intersect the zero trajectory curve, the control effort 
changes, according to (2.43) to u = +N, and follows the zero trajectory curve 
into the origin. 

2. Limit Cycles 

The control effort, u = +N will not only drive the states to the origin, it 
will, in fact, become confused there. It is a point of indecision and chatter or 
limit cycle motion ensues as in Figure 2.12. The magnitude of the limit cycle 
depends on the sampling rate or time delay for the control effort. If the sample 
rate is low, the states will penetrate far into the opposite control region before 
the control effort can change, and the limit cycle will swing widely around the 
origin. If the sample rate is high then the states will not travel far away from 
the origin before the control effort corrects the direction of travel. 

Since a heavy chatter mode is not desirable for many real systems, 
control logic for reducing or removing the chatter mode may be developed, 
however, it will not be included in this paper. 
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Figure 2.12 Limit Cycle on Second Order Solution Trajectory 
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3. Simulation of the Switching Law Controller 

To test the switching law (2.43), we simulate the system of (2.29) 
using a maximum control effort of N = 1 and initial conditions Xi(0)>0, 
X2(0) > 0. The output of the simulation, shown in Figure 2.13, demonstrates the 
control effort driving the states first to the zero trajectory curve, then along the 
zero trajectory curve to the origin. Once at the origin, the control effort goes 
into limit cycle as shown in Figure 2.14. 
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Figure 2.13 Simulation Using Second Order Switching Law 




Figure 2.14 Control Effort for Second Order Switching Law 
Solution 
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III. APPLICATION OF A SECOND ORDER CONTROLLER 



A. MISSILE/TARGET INTERCEPT MODEL 

An application for Minimum Time Optimal Control is an anti-missile defense 
system, where the incoming target has a speed advantage over the defensive 
missile. In this case the missile control system must operate in saturation mode. 
The bang-bang controller, where control effort is either maximum-positive or 
maximum-negative, has a faster response capability than the standard 
Proportional Navigation guidance system. A simple model of a missile to target 
engagement can be described by Figure 3.1 where a is the line of sight (LOS) 
angle from missile to target, Ym is the angle of the missile velocity, yt is the angle 
of the target velocity, and all angles are relative to an inertial reference. 




Figure 3.1 Missile/Target Geometry 
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It is assumed that the target is on the final leg of it’s flight and is now on a 
straight, non-maneuvering trajectory. The control is to drive the line of sight 
rate (a) and it’s derivative (a) to zero in minimum time. 

For our models we will be using only two dimensional scenarios. The 
system dynamics use second order models for each dimension. The target 
dynamics are non-maneuvering so that the acceleration, at, is zero. The missile 
acceleration, a^, is assumed perpendicular to the velocity vector Yni- The 
system dynamics are shown in the signal flow graph in Figure 3.2. 
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The LOS angle, a, and it’s derivatives, a and o, may be derived analytically 
from the available states. We define the geometry of the system as in Figure 
3.3. 




Figure 3.3 Geometry for Angle Definitions 
with 



Ri = Range of the Target. 

Rm = Range of the Missile. 

Vi = Velocity of the Target. 

Vm = Velocity of the Missile, 
at = Acceleration of the Target. 
Um = Acceleration of the Missile. 
Vc = Closing Velocity. 

so that 



o = tan 



r R ^ 

^ ^^imx y 



(3.1) 
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D .V -R .y 
• _ *^imx ’imy *^tmy ^imx 



(3.2) 



and 



.. _ ^tmx ' ^tmy ^^tmx ' ^tmy ^tmy'^tmx ^ ^ ^ tmx ' ^tmy ^tmy ' ^tmx ) .^ 

j^2 j^3 ^ ) 

where 

V, = -R (3.4) 

In an actual application Kalman or Luenberger Observers may be used to 
obtain estimates of 6 and d. 



B. APPLICATION OF THE SWITCHING LAW 

The switching law from (2.43) is applied to our scenario where d and d 
form the state space of this simulation, and is implemented as 



u = N • sign 



r 



a + 



aja 



(3.5) 



The sign convention of (2.43) was changed to conform to the geometry of 
scenario. 

We set the simulation with the initial conditions 



our 



Rmx(O) = 0 ft 
Vn,x(0) = 2000 ft/sec 
amx(O) = 0 ft/sec2 
Rtx(O) = 10000 ft 
Vtx(O) = 2000 ft/sec 
aix(O) = 0 ft/sec2 
^finai — 2.25 sec 



Rmy(0) = 0ft 
Vn,y(0) = 0 ft/sec 
amy(O) = 0 ft/Sec2 
Rty(O) = 1000 ft 
V,y(0) = 0 ft/sec 
aty(O) = 0 ft/sec2 
dt = 0.01 sec 
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Running the simulation we find that the controller does drive a and d to 
zero, and maintains them about the origin until intercept is reached (Figure 3.4). 
Once the states reach the origin, the system controller goes into chatter mode, 
see Figure 3.5, and the trajectory chatters back and forth about the desired 
intercept path. 
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Figure 3.4 Missile/Target Intercept Simulation Using Second 
Order Controller 
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Figure 3.5 Control Effort for Second Order Missile/Target 
Simulation 
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C. COMPARISON WITH PROPORTIONAL NAVIGATION 

CONTROLLER 

The system dynamics for the proportional navigation controller simulation is 
basically the same as that in the previous section, except that the control effort , 
u, is proportional to a, which is obtained with an estimator, shown in Figure 3.6, 
where P = d. We limit the acceleration so that u is bounded by ±N. With only 
the change in the controller we ran the simulation for the same initial conditions 
with the results shown in Figures 3.7 and 3.8. 

This Proportional Navigation system is effective at intercepting the target 
only when the missile/target geometry and kinematics are sufficient that the 
missile has time to maneuver. The delay in saturation of the control effort 
caused the missile to turn too slowly so that the Proportional Navigation 
Controlled system was unable to maneuver quickly enough to intercept the 
faster, although non-maneuvering, target. 
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Figure 3.6 System Dynamics of Proportional Navigation 
Missile/Target Model 
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Figure 3.7 Proportional Navigation Missile/Target Simulation 
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Figure 3.8 Control Effort for Proportional Navigation 
Missile/Target Simulation 
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IV. THIRD ORDER CONTROLLER 



A. FORWARD TIME SYSTEM 

The second order controller is effective at intercepting a target that has 
speed advantage. Previously we have controlled the system by driving <j and 
d to zero, and thereby maintaining a constant LOS angle, a, until impact. But 
what about controlling the LOS angle itself? There are situations in which we 
would like to be able to attack a target from a particular angle, or perhaps have 
multiple missiles, each with it's own pre-defined attack angle. 

Consider a point defence system in which the missile, launched from some 
point away from the target's final trajectory, would try to position itself in 
minimum time onto a "head-on" collision course with the target, as shown in 
Figure 4.1. In such a situation we would use a third order minimum time 
controller to drive a, d, and d to zero. 

1. System Definition 

Our third order example is 





‘0 
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0" 




‘O' 


X = 


0 


0 
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x-t- 


0 
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0 


0 _ 




1 



Minimizing the Hamiltonian and solving the system we find 

H = 1 -I- PjX2 + P2X3 -f- P3U 
and 

u = -N-sign(p3) 



(4.1) 



(4.2) 

(4.3) 
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Attacking Missile 




Figure 4.1 Application of Third Order Minimum Time Controller 
where 



so that 




0 0 

0 -1 

0 0 



0 

0 

-1 



Pi =Ci 

P2 “ C|t " 1 “ C2 

P3 =iCit^-C2t + C3. 



( 4 . 4 ) 



( 4 . 5 ) 

( 4 . 6 ) 

( 4 . 7 ) 
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From the adjoint solution the control may switch sign no more than 
twice. Again tracing this problem in negative time we may follow the zero 
trajectory curves out from the origin with control efforts of ±N, These curves 
are now in three dimensional space residing in their own plane which intersects 
the origin. Intersecting these zero trajectory curves are an infinite number of 
curves making a surface, and leading off from this surface the infinite number of 
trajectories lead to the initial conditions. Therefore in forward time we may 
start with an initial condition such that u = +N will drive the system to intersect 
with the surface at tg^i as shown in Figure 4.2. Switching the control effort to 
u = -N will drive the system along the surface to intersect with the zero 
trajectory curve at tsw2 where the control effort switches again. Finally u = +N 
will drive the system to the origin. 




Figure 4.2 Three Dimensional Switching Curves 
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2. Third Order Switching Curves 

We have determined from the second order system that solving for the 
control switching times is not as widely applicable as defining the control as a 
function of the states; so we now move on to defining the third order switching 
curves. 

Starting with the state equations (4.1) we discretize and expand the 
equations so that 
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A3 is the zero matrix so it and all higher terms of A are dropped. 
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Combining these we obtain 



(4.8) 



(4.9) 
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It^l 
2 ^ 






x(t) = 


0 
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x(0) + 






0 
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t 



or in scalar equations 

X,(t) = Xi(0) + tX2(0) + it2x3(0) + lt^u(0) 
X2(t) = X2(0)+ 1X3(0) + ^t^u(O) 

X 3 (t) = X3(0)+tu(0). 



(4.10) 



(4.11) 

(4.12) 

(4.13) 



B. NEGATIVE TIME SYSTEM 

The third order system, being piecewise continuous, is easily broken into 
several simple boundary value problems. In order to solve the systems of 
curves we must determine some boundary values for the intersections of these 
families of curves. We start at the origin and run the system in negative time to 
determine our other boundary conditions. 

In negative time we have 



Discretizing we find 
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(4.14) 



(4.15) 
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(4.16) 



so that 
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Jo 



Bdx = 



-it" 

-t 





’l 








x(t) = 
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-t 


x(0) + 
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0 


1 




-t 



u(0), 



(4.17) 



or in scalar equations 

x,(t) = x,(0)-lX2(0) + it%(0)-lt\(0) (4.18) 

X2(t) = X2(0)-IX3(0) + 4t^u(0) (4.19) 

X 3 (t) = X3(0)-tu(0). (4.20) 

1. Solving for Negative Time Boundary Points 

To develop a complete solution we solve the equations for several 
different points along the zero trajectory curve. Setting u = -N and traveling out 
along the zero trajectory curve for 1 second we find 



x,(l) = -iuc’=-l(-N) = iN 


(4.21) 


X 2 (l) = 4ut^=4(-N) = -4N 


(4.22) 


X 

U) 

II 

1 

c 

II 

1 

1 

II 


(4.23) 



Similarly we run the system in negative time for 2 and 3 seconds to 
obtain other boundary points as shown in Figure 4.3, and listed in TABLE 1. 
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NEGATIVE TIME BOUNDARY CONDITIONS, u = -N 



u = -N 


t = 2 


t = 3 


Xl(t) 


IN 


|N 


X2(t) 


-2N 


-|N 


X3(t) 


2N 


3N 



TABLE 1. 



The control effort may also be u = +N, traveling away from the origin 
on the other zero trajectory curve giving the boundary conditions in TABLE 2. 
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NEGATIVE TIME BOUNDARY CONDITIONS, u = +N 



u = +N 


t= 1 


t = 2 


t = 3 


Xl(t) 


-iN 


jN 


|N 


X2(t) 


IN 


-2N 


-fN 


X3(0 


-N 


2N 


3N 



TABLE 2. 



We may now solve the equations for the family of curves that intersect 
the zero trajectory curves. Specifically, we solve for the equation of the one 
curve that travels from some point xi(0), X 2 ( 0 ), X 3 ( 0 ), to the point xi(t) = -^N, 
X 2 O) = "jN, X 3 (t) = N. Since the control effort on the zero trajectory curve for 
this intersection point is u = -N, the control effort for the curve we are solving 
for must be u = +N, and the forward time equations may be written as 

x,(t) = lN = x,(0)+tX2(0) + il%(0) + it’N (4.24) 

X2(t) = -|N = X2(0)+tX3(0)+it’N (4.25) 

X3(t) = N = X3(0) + tN (4.26) 

Solving (4.26) for t we find 



t = 



X3(0) 

N 



+ 1 . 



(4.27) 



Substituting (4.27) into (4.24) and (4.25), and after a bit of algebra we obtain 



-N = x2(0)- 



X3^(0) 

2N 



(4.28) 



0 = Xi(0) + X2(0)- 



X2(0)X3(0) Xj^O) , X33(0) 



+ ■ 



(4.29) 



N 2N 3N^ 

Using the other boundary values from TABLE 1 and TABLE 2 as 
well, we may generate a family of equations representing both sides of the zero 
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trajectory curves, and different points of intersection along the curves (see 
TABLE 3). The control effort, u, in TABLE 3. is the control effort for the zero 
trajectory curve that is intercepted. The time, t, is the time out from the origin 
to the intercept point for the negative time system. 

To combine all the equations from TABLE 3 into one solution we 

define 



( 



X-, • 



w = sign 



Xo +■ 



V 



2N 



and 



f = Xo + w 



so that our family of curves is defined as 



2N 



(4.30) 



(4.31) 



0 = xi(0) + 



X3 (0) , X2(0)-X3(0) 



3N^ 



+ w 



N 




(4.32) 



Equation (4.30) determines which signs are to be applied based on which 
direction the zero trajectory curve is on; Equation (4.31) adjusts the magnitudes 
of the equation depending on the distance of the intersection of the zero 
trajectory curve from the origin. Each curve intersecting the zero trajectory 
curve will have a different value of the function f, and f will remain constant all 
along that curve. 
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FAMILY OF SOLUTION CURVES 



u = -N 


N = x2(0) 

2N 


t = 1 


* ^ N 2N 3N^ 


u = -N 


4N = x2(0) 

2N 


t = 2 


0 = x,(0) + 2x,(0)- 

^ ^ N N 3N^ 
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t= 1 


0 = x,(0) + x,(0)+ V( 0 ) 

* ^ N 2N 3N^ 
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II 
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4N = x2(0)+^^^^^^ 
2N 


t = 2 


0 = x,(0) + 2x,(0)+ ^ X 3 ^( 0 ) ^ X33(0) 

’ ^ N N 3N^ 


u = +N 


9N=x2(0)+^^ 

2N 


t = 3 


0 = x.(0).3x3(0).’'^(0)’^3(0),3x3^0)^X33(0) 

* ^ N 2N 3N^ 



TABLE 3. 
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C. THIRD ORDER SWITCHING LAW 



Using (4.30), (4.31) and (4.32) we can break up space into areas of 
opposite control effort. The control effort is defined by which of these volumes 
contain the states. Therefore, the third order switching law is defined as 



u = - N - sign 



xi(0) 



. X3^0) 

3N^ 



+ 



w 



X2(0)-X3(0) ^ ^ 

N 




(4.33) 



D. THIRD ORDER CONTROLLER SIMULATION 



A simulation of this third order model shows that the third order switching 
law, (4.33), drives the states to the origin with only 2 changes is the control 
effort. In application some means of shutting off the control effort will be 
required for the system to avoid entering chatter mode upon reaching the origin 
(see Figure 4.4). 

An examination of the functions (4.30), (4.31), and (4.32) in Figure 4.5 
shows that f is a smoothly increasing curve until the states intercept the 
switching surface. Once the states are following the surface, f remains a 
constant value. When the states reach the zero trajectory curve, f becomes 
approximately zero. 
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Figure 4.4 Simulation of a Third Order Minimum Time 
Controller From a Fixed Initial Condition 
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Figure 4.5 Control Laws For Third Order Solution 
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V. MISSILE ADJOINTS 



Given a time varying system 

X = A(t)x + B(t)u 
y = C(t)x. 

it can be shown that the impulse response of the adjoint system 



(5.1) 

(5.2) 



p = A(tf-t)'^p + C{tf-t)^r 
y = B{tf-t)’^p 



(5.3) 



(5.4) 



is the response of the original system, at time t, to an impulse applied at time 
tf-t before t. [2] For a time invariant system the transfer function for the 
original system is identical to the transfer function of the adjoint system, i.e., 
they are self-adjoint. 

A. CONSTRUCTION OF AN ADJOINT 

Given the mathematics of the adjoint method, we now need to define some 
rules to create and understand the adjoint system with real systems and block 
diagrams. [3] 

1. Rule 1: Convert All System Inputs to Impulses 

In constructing the adjoint it is necessary that all system inputs be 
impulsive. Since this may not be the case with the block diagram, all system 
inputs and initial conditions must be converted to impulsive inputs via block 
manipulations and extensive use of integrators.. Figure 5.1 shows that step 
inputs and initial conditions are equivalent to the output of impulse driven 
integrators. 
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2. Rule 2: Replace t With tf - t in the Arguments of All Time 

Varying Coefficients 

In many linear systems it is possible to express a gain as a function of 
time. The adjoint system operates in negative time and Figure 5.2 shows the 
conversion of functions of time into the adjoint domain. 






Oricinal Svstem 


Adjoint Svstem 


Time 


K(t) = at+ b 


K(t-tf) =a(t-tf) + b 


Varying 

Gain 


K(t) = -. \ 

a(t-tf ) + b 


K(t-t,) = — !— 

at +b 



Figure 5.2 Convert Functions of Time to the Adjoint Domain 
3. Rule 3: Reverse the Direction of all Signal Flow 

The direction of all signal flow must be reversed, redefining nodes as 
summing junctions and visa versa. Notice that all system outputs become 
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inputs and all system inputs become outputs. This last rule allows for the simple 
graphical creation of the adjoint system by first drawing the block, or flow, 
diagram of the original system, then redrawing it with the arrows reversed. 
Figure 5.3 shows some examples f converting nodes to summing junctions and 
visa versa. 



Original System 



Adjoint System 



G(s) 

-> 



h b(t) 

>► 






G(s) 



Yb(t-tf ) 

— < 



Figure 5.3 Converting Nodes to Summing Junctions 



B. DEVELOPMENT OF A SECOND ORDER ADJOINT SYSTEM 
This adjoint formulation lends itself well to analyzing some optimization 
problems, those where tf is free and the terminal state is constrained to a point, 
curve or surface. From Pontryagin the control is a function of the homogeneous 
adjoint (see (4.2) and (4.3)). The missile adjoint form gives also a particular 
solution of the adjoint (i.e. system impulse response). 

The adjoint allows one to generate optimum solutions (i.e. switching 
surfaces) for all possible initial conditions. The forcing impulses passed through 
an integrator gives our saturation type of optimal control (±N). 

As an example for the application of the method of adjoints we will apply 
our adjoint rules to our second order controller. The system described by (2.9) 
has a time dependent control input where 
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(5.5) 



u(t) = 




t<ts 

t>ts 



and, from our second order example, 



ts = 



/x(0) l/"x(0)^^ 



N 



+ - 

2 



y N J 



+ 



x(0) 

N 



(5.6) 



1. Apply Rule 1 to the Second Order Example 

For this example we will define N = 1. We first draw out the original 
signal flow diagram, showing all the inputs and initial conditions. In order to 
apply the impulsive inputs we expand the second order system to the third 
order state equations 





’0 
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o’ 




'o’ 


X = 
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x + 
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0 




1 



y = [l 0 0]x. 

Converting all the inputs to impulsive inputs and 
flow diagram in Figure 5.4. 



u (5.7) 

(5.8) 

initial conditions we get the 
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Figure 5.4 Application of Rule 1 on Second Order Example 
2. Apply Rule 2 to the Second Order Example 

We replace t with t - tf in the arguments of all time varying coefficients 
to obtain the flow diagram in Figure 5.5. 
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3. Apply Rule 3 to the Second Order Example 

Finally we reverse the direction of the signal flow, redefining nodes as 
summing junctions and visa versa, thereby changing the system inputs to 
outputs and the system outputs to inputs as shown in Figure 5.6. 




The mathematical definition of the Adjoint System is 

P=a'^P + C'^u (5.9) 

y = B’^p. (5.10) 

Using the A, B, and C matrix from (5.7) and (5.8) we obtain 
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_0 


1 


0 _ 


.P3. 




0 _ 



y = [0 0 l]p 

which corresponds to Figure 5.6. 



(5.11) 

(5.12) 



C. SIMULATION OF THE SECOND ORDER ADJOINT 
1. Forward Time Second Order Simulation 

In the forward time second order simulation we set N = 1, and the 
initial conditions x(0) = 2 and x(0) = 0. Using impulses through integrators we 
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apply the initial control effort at t = 0. A second impulse, opposite in sign and 
twice in magnitude, is applied at the switching time, as defined in (5.6), driving 
the states through the origin. 

The simulation length is 4 seconds with a sample step of .05 seconds. 
The results are shown in Figure 5.7. Starting at the initial conditions the system 
is driven through the origin with a minimum miss distance of 0.003344 at a time 
of 2.83 sec. 

2. Adjoint Solution for Second Order System 

In the adjoint domain the system will travel in negative time starting at 
the origin and traveling outward to the initial conditions for the forward time 
system. Because this is a time invariant system the trajectory for the adjoint 
solution will be the same as the forward time system but in the opposite 
direction. From our second order example, the optimum switching in negative 
time from the terminal state at the origin is 



The trajectory constraint yields 




(5.13) 







X2(ts)|x2(ts)| 

2N 



(5.14) 



A reverse impulse of twice the magnitude is applied and the trajectory 
proceeds out to the desired initial conditions given. The negative time impulse 
response is thus prescribed by using the switching times 



tf-t,=lVNxi(0) + lx2(0)' 



(5.15) 
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t,w-to = ■^VN>‘,( 0 ) + ix2(0)^ 



(5.16) 



The parameters for the adjoint solution are the same as for the 
forward time, but the initial values for the states are at the origin. The output 
of the adjoint system, according to (5.12), is ps, so the phase plot of Figure (5.8) 
plots -p 2 and p 3 . The adjoint solution traces almost the same path as the 
forward time solution, switching at tj = 1.582 sec., and missing the point (2,1) by 
0.002733 at tf = 2.83 seconds. By changing the initial sign of the control effort, 
and adjusting the switching time, we could drive the adjoint system to any 
desired point in the phase plane, and therefore know the optimal solution for the 
forward time system. 



49 




Figure 5.7 Simulation of Forward Time System 




Figure 5.8 Simulation of the Adjoint Solution 
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D. MISSILE SIMULATION WITH ADJOINTS 



The method of adjoints can be extremely useful when dealing with time 
varying systems. By using the adjoint solution we are able to see how the 
forward time system behaves for all times. 

1. Missile/Target Model 

We will simplify our scenario of missile/target engagement in order to 
present an uncluttered example of the method. The target will have zero x- 
velocity and constant y-velocity, and will start on the x-axis at a distance, R, 
from the origin (see Figure 5.9). The missile will start at the origin with a 
constant x-velocity by zero initial y-velocity. 




Figure 5.9 Geometry of Missile/Target Adjoint Solution 

Since R is large and yi - y„i is small we use the small angle 
approximation and define 



o = tan 



-1 



I R J 



- Yl - Vr 

R 



The closing velocity is constant so we note that 

R = V,(tf-t) 



(5.17) 



(5.18) 
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where tf is the length of the simulation. 

2. Signal Flow Diagram 

Because the x-velocities are constant we need only to model the y- 
dimension. We develop the flow diagram in Figure 5.10 from the geometry of 
Figure 5.9. We use proportional navigation control and impulsive inputs are 
used for initial conditions. 




The estimator used was developed from the mechanics and dynamics 
of a seeker head system, however, it can be shown to be equivalent to a 
Kalman Estimator or Luenberger Observer. This estimator in Figure 5.1 1 has a 
time constant of 0.1 seconds. The primary interest in this model is to determine 
the miss distance at the final time, t = tf, so the output of the system is ymiss- 
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Figure 5.11 Signal Flow Diagram of the Estimator 
3. Forward Time Simulation 
We define our states to be 
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and the system state equations are then 
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(5.20) 



0 0 0 0 0 



0 



y = [-l 0 0 0 1 0 ]x. 



(5.21) 
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Since the x-velocities are constant we know approximately when the 
point of closest approach will occur, and therefore about how long to run the 
simulation. The LOS angle, a, is a function of t and tf, so that different values of 
tf will result in changes to the control effort and system response. Because this 
state matrix. A, is time dependent , it must be redefined at each time step of the 
simulation. We define our proportionality constant, n, to 4, the closing velocity, 
Vc, to 5000 ft/sec., so that with tf = 4 sec. the initial range is 20,000 ft. The 
output of the simulation is shown in Figure 5.12. 

To study the effect of different values of tf, or to locate the optimal 
value, we may have to run the forward time simulation many times, which 
could be a very tedious process, especially since we are only interested in the 
final value of t = tf. 

4. Adjoint Solution 

An alternative to running the forward time system over and over again 
would be to run the adjoint solution once to generate the final values of the 
family of forward time solutions. This time we generate the adjoint system 
using matrix algebra instead of the graphical method. The solutions to the 
adjoint system are, from (5.9) and (5.10) 
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(5.22) 



y = [0 0 0 0 0 l]p. (5.23) 

This system generates a curve whose value at any time, tg, is the final 
value of the forward time system where tf = tg. Figure 5.13 shows four curves 
from the forward time solution corresponding to tf = 1,2,3, and 4 sec. Each 
curve end exactly on the adjoint solution curve at the corresponding adjoint 
time. 
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Figure 5.12 Forward Time Simulation 




Figure 5.13 Adjoint Solution Curve 
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VI. CONCLUSIONS AND COMMENTS 



We have developed the second order minimum time optimal controller and 
applied it to the fast reaction missile defense problem with both the closed form 
analytic solution, and the open form switching curves solution. In comparing 
our open form solution with a standard proportional navigation controller 
solution we demonstrate an increased maneuverability enabling our missiles to 
defend against faster, more capable threats. The accuracy at intercept is a 
function of the control logic used to shut off the control effort when desired 
conditions are met. This is a subject that should be explored in future projects. 

We introduced a third order minimum time controller which promises to not 
only improve reaction time and maneuverability, but also presents us with the 
ability to control and define a desired attack angle for an improved destructive 
potential. A model for a practical third order controller should be developed 
and evaluated. 

Having introduced the method of adjoints and shown some of its functions, 
we hope to stimulate some more practical applications of this technique. We 
are excited at the possibilities of using the method of adjoints in determining 
optimal, closed form solutions to forward time problems at speeds fast enough 
for practical implementation. 
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APPENDIX A. PROGRAM CODE 



All of the simulations for this project were run on both IBM-ATi class and 
Macintosh IP computers using the matrix manipulation language MATLAB^. 
For IBM-AT compatibles MATLAB, version 3.5f was used, and on the 
Macintosh II computers version MAC II-MATLAB, version 1.1b. This 
appendix contains the source code for all of the simulations and functions 
written in support of this project. 

Only a limited background in programming is required for understanding 
these files. While MATLAB is similar to FORTRAN, MATLAB's control 
structures are much less complex, and with matrix manipulation built into the 
system, vector definition and storage are greatly simplified. Comments are 
started by the percent sign (%) and continue to the end of the line. Ellipsis (...) 
at the end of a line indicate the continuance of the logical line onto the next line 
of code. 

Each file will begin on a new page to assist those who are interested in 
examining or reproducing the code. Although an analysis of these files is not 
necessary to understand this report, readers are encouraged to examine them 
closely for further information. 

The code is presented in the order of usage in the main text with all 
supporting functions grouped with the main program of interest. 



1 IBM and IBM-AT are registered trademarks of IBM. 

^ MAC II is a registered trademark of APPLE. 

^ MATLAB is a registered trademark of The MathWorks, Inc. [4] 
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1. BB2NDST.M 



% BB2NDST.M 11 Mar. 1991 

% BB2NDST.M is a simulation of the 2nd order bang-bang controller with an analytic 
% solution for the Switching Time of the control effort. The signs of the control effort 
% must be matched with the initial conditions to have convergence. 

% written by Colin R. Cooper 

% Define the State Equations. 

A = [0 1;0 0]; 

B = [0 1]'; 

C = [10]; 
xl0 = 2; 
x20 = 1 ; 

N = 1; % Maximum control effort 
aN = abs(N); 

Tf = 4.17; 
dt = .005; 

[Phi,Del] = c2d(A,B,dt); 

% Create the storage vectors, 
kmax = Tf/dt -i- 1; 

X = zeros(2,kmax); 
y = zeros(l,kmax); 
time = zeros(l,kmax); 

x(:,l) = [xl0;x20]; % Set initial conditions for x. 

% Define the Switching Time for the control effort, 
tsf = abs(x20)/aN + sqrt(abs(xlO -t- x20*abs(x20)/(2*aN))/aN); 

% Begin simulation loop. 



% xl initial condition. 

% x2 initial condition. 

% Length of simulation. 

% Time increment for simulation. 
% Discretized System. 
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for i = l:kmax - 1 
if time(i)<tsf 
u = -N; 

else 

u = N; 
end 

x(:,i+l) = Phi*x(:,i) + Del*u; 
y(l,i+l) = C*x(:.i+1); 

time(i+l) = time(i)+dt; 
end 

% Plot the output of the simulation. 
clg,plot(x(l,:),x(2,:),'-w'),grid 
title(’Phase plot - Switching Time’) 
xlabel('Xl’),ylabel(’X2’) 



2. BB2NDSL.M 



% BB2NDSL.M 11 Mar. 1991 

% is a simulation of the 2nd order bang-bang controller using a Switching Law for the 
% control effort. 

% written by Colin R. Cooper 



A = [0 1;0 0];% 
B = [0 1]’; 

C = [10]; 


Define the State Equations. 


xlO = 2; 


% xl initial condition. 


x20 = 1; 

N = 1; % Maximum control effort. 


% x2 initial condition. 


Tf = 4.5; 
dt = .01; 

[Phi,Del] = c2d(A,B,dt); 


% Length of simulation. 

% Time increment for simulation. 
% Discretize the System. 


% Create storage vectors, 
kmax = Tf/dt + 1; 

X = zeros(2,kmax); 
y = zeros(l,kmax); 
u = zeros(l,kmax); 
x(:,l) = [xl0;x20]; 


% Set initial conditions for x. 


% Begin simulation loop, 
for i = 1 :kmax - 1 





u(i) = -N*sign(x(l,i) + .5*x(2,i)*abs(x(2,i))/N); 
x(:,i-i-l) = Phi*x(;,i) + Del*u(i); 

(l,i-t-l) = C*x(:,i+1); 
end 
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% Plot the output of the simulation. 
clg,plot(x( 1 ,:),x(2,:)),grid 
title('Phase plot - Switching Law') 
xlabel(’Xr),ylabel('X2') 
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3. SIM2SL.M 



% SIM2SL.M 2nd Order Switching Laws Control 08 Apr. 1991 

% Simulation of the missile/target simulation. 

% 

% -> Uses analytic values for sigma-dot,sigma-ddot for calcultion of the control effort 
% using the 2nd order switching curves. 

% -> The Target makes no evasive maneuvers. 

% -> Calculates the Quantization Error based on the average velocities for the crossover 
% endpoints. 

% -> Calls INTERP.M function which takes the states at the crossover endpoints, creates 
% 100 point straight lines to connect the points, and determine the minimum miss 

% distance. 

% -> Allows delay time for missile control. 

% -> Target is now on level flight with Beta = 0. 

% written by Colin R. Cooper 



% Define states. 



N = 1000; 


% Maximum control effort. 


Am = [0 1 0 0;0 0 0 0;0 0 0 1;0 0 0 0]; 
Bm = [0 0;10;0 0;0 1]; 


% Missile State Equations. 


At = [0 1 0 0;0 0 0 0;0 0 0 1;0 0 0 0]; 
Bt = [0 0;l 0;0 0;0 1]; 


% Target State Equations. 


Tf=2.25; 


% Total time of simulation. 


dt = .01; 


% Sample step size. 


% Descretize the states. 


[Phim,Delm] = c2d(Am,Bm,dt); 


% Discrete Missile System 


[Phit,Delt] = c2d(At,Bt,dt); 


% Discrete Target System 
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% Define storage vectors, 
kmax = Tf/dt + 1; 
xm = zeros(4,kmax); 
xt = zeros(4,kmax); 



% xm = [ y yd X xd ]’ 
% xt = [ y yd X xd ]' 



time = zeros(l,kmax); 
sig = zeros(3,kmax); 
am = zeros(l,kmax); 
um = zeros(2,kmax); 

Rtm = zeros(l,kmax); 

xm(:,l) = [0;0;0;2000]; % Initial Conditions for Missile. 

xt(:,l) = [ 1 000; 0; 10000; -3 000]; % Initial Conditions for Target. 

Rtm(l) = sqrt((xt(l,l)-xm(l,l))^2 + (xt(3,l)-xm(3,l))^2); % First value for Range, 
acm = 0.0; % Initial acceleration for the missile, 

act = 0.0; % Initial acceleration for the Target. 

% Begin the Simulation Loop, 
for i = l:kmax-l 
% Define angles. 

beta = atan2(xt(2,i), xm(4,i)); % Velocity angle for the Target, 

gam = atan2(xm(2,i), xm(4,i)); % Velocity angle for the Missile. 

sig(l,i) = atan2(xt(l,i)-xm(l,i), xt(3,i)-xm(3,i)); % LOS angle. 

sig(2,i) = ((xt(3,i)-xm(3,i))*(xt(2,i)-xm(2,i)) - (xt(l,i)-xm(l,i))... 
*(xt(4,i)-xm(4,i)))/Rtm(i)'^2; 

sig(3,i) = ((xt(3,i)-xm(3,i))*(act*cos(beta)-acm*cos(gam))-... 
(xt(l,i)-xm(l,i))*(act*sin(beta)+acm*sin(gam))+... 
2*(xt(4,i)-xm(4,i))*(xt(2,i)-xm(2,i)))/Rtm(i)'^2; 
sig(3,i) = sig(3,i)+2*(xt(4,i)-xm(4,i)+xt(2,i)-xm(2,i))*sig(2,i)/Rtm(i); 

% Acceleration = Max value perpendicular to gamma. 
am(i) = N*sign(sig(2,i) + sig(3,i)*abs(sig(3,i))/(2*N)); 
um(:,i) = [am(i)*cos(gam); -am(i)*sin(gam)]; 
xm(:,i+l) = Phim*xm(:,i) + Delm*um(:,i); 
xt(:,i+l) = Phit*xt(:,i) + Delt*[0;0]; 
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time(i+l) = time(i) + dt; 

Rtm(i+1) = sqrt((xt(l,i+l)-xm(l,i+l))^2 + (xt(3,i+l)-xm(3,i+l))'^2); 
acm = am(i); 



end 



% Evaluation of Miss distance, 
iflag = 0; 

rm = find(Rtm==min(Rtm)); 
if rm == kmax 



% Index of the minimum range value. 

% Determine whether the crossover occurs 
% before or after the min range value. 



iflag = 1; 



It = rm - 1 ; 

elseif Rtm(rm+ 1 )<Rtm(rm- 1 ) 

It =rm; 
else 

It = rm - 1; 
end 

% Average Velocities in Intercept Area. 

Vm = .5*(sqrt(xm(2,lt)^2+xm(4,It)^2) + sqrt(xm(2,It+l)^2+xm(4,It+l)^2)); 

Vt = .5*(sqrt(xt(2,It)^2+xt(4,It)^2) + sqrt(xt(2,It+l)^2+xt(4,It+l)^2)); 

QE = .5*dt*(Vm + Vt); % Quantization Error, 

if iflag == 0 

r = interp(xm(:,It:It+l),xt(:,It:It+l)); % Interpolated minimum miss distance, 
rstr = ['Inter Miss = ' num2str(r) ' ft']; 
else 

rstr = ['Crossover never reached']; 
end 

% Plot the output of the simulation with results. 

plot(xm(3,:),xm(l,:),'-w',xt(3,:),xt(l,:),'-w') 

text(.15,.85, ['Intercept Time = ' num2str(time(rm)) ' sec'],'sc'); 

text(.15,.81,['Min. Miss = ' num2str(Rtm(rm)) ' ft'], 'sc'); 

text(.15,.77,['Quan. Err = ' num2str(QE) ' ft'],'sc'); 

text(.15,.73,rstr,'sc'); 
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text(. 15,.69,['N = ',num2str(N)],'sc'); 
title('Control; 2nd Order Switching Laws (sig, sigd)’) 
xlabeK'Direction 1 (ft)'),ylabel(’Direction 2 (ft)') 
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4. SIMPRNV.M 



% SIMPRNV.M simulation of the missile/target simulation. 
% Proportional Navigation Guidance. 

% written by Colin R. Cooper 



% Define states. 

As = [0 1; -64-16]; 

Bs = [0; 64]; 

Am = [0 1 0 0;0 0 0 0;0 0 0 1;0 0 0 0]; 
Bm = [0 0;1 0;0 0;0 1]; 

At = [0 1 0 0;0000;000 1;0000]; 
Bt = [0 0;1 0;0 0; 0 1]; 

Tf = 2.25; 
dt = .01; 
maxac = 1000.0; 



% Estimator for sigma-dot. 

% Missile State Equations. 

% Target State Equations. 

% Simulation Time. 

% Sample step size. 



% Descretize the states. 

[Phis, Dels] = c2d(As,Bs,dt); 
[Phim,Delm] = c2d(Am,Bm,dt); 
[Phit,Delt] = c2d(At,Bt,dt); 

% Create strage vectors, 
kmax = Tf/dt -i- 1; 
xm = zeros(4,kmax); 
xt = zeros(4,kmax); 
um=zeros(2,kmax); 
time = zeros(l,kmax); 

Rtm = zeros(l,kmax); 
b = zeros(2,kmax); 



% Estimator System. 
% Missile System. 

% Target System. 



% xm = [ y yd X xd ]' 
% xt = [ y yd X xd ]' 
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xm(:,l) = [0;0;0;2000]; % Missile Initial Conditions. 

xt(:,l) = [1000;0;10000;-3000]; % Target Initial Conditions. 

Rtm(l) = sqrt((xt(l,l)-xm(l,l))^2 + (xt(3,l)-xm(3,l))'^2); % First value for Range. 



% Begin simulation loop, 
for i = l;kmax-l 



Vm = sqrt(xm(2,i)'^2 + xm(4,i)^2); 

ub(i) = atan2(xt(l,i)-xm(l,i), xt(3,i)-xm(3,i)); 



b(:,i+l) = Phis*b(:,i) + Dels*ub(i); 
gam = atan2(xm(2,i), xm(4,i)); 
if Vm*b(2,i)<=maxac 



% Apply Thrust limitation. 



am = Vm*b(2,i); 

else 

am = maxac; 
end 

um(:,i) = [am*cos(gam); -am*sin(gam)]; % Acceleration is perpendicular to gamma. 
xm(:,i+l) = Phim*xm(:,i) + Delm*um(:,i); 
xt(:,i+l) = Phit*xt(:,i) + Delt*[0;0]; 
time(i+l) = time(i) + dt; 

Rtm(i+1) = sqrt((xt(l,i+l)-xm(l,i+l))^2 + (xt(3,i+l)-xm(3,i+l))^2); 



end 

rm = find(Rtm==min(Rtm)); 
if Rtm(rm+ 1 )<Rtm(rm- 1 ) 



% Index of minimum range value. 

% Determin whether crossover occurs 
% before or after the min range value. 



It = rm; 

else 



It = rm - 1; 
end 

r = interp(xm(:,It;It+l),xt(:,It;It+l)); 
Vm = sqrt(xm(2,It)^2 + xm(4,It)^2); 
Vt = sqn(xt(2,It)^2 + xt(4,It)^2); 

QE = .5*dt*sqrt(Vm'^2 + Vt'^2); 



% Quantization Error. 



% Obtain the Interpolated min miss distance 



% Plot the output of the simulation. 
axis([0 10000 0 1400]); 



% Same scale as the Optimal Control Sim. 
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plot(xm(3,:),xm(l,:),'-w',xt(3,:),xt(l,;),'-w') 
text(. 15, .85, [’Intercept Time = ’ num2str(time(rm)) ' sec'],' sc') 
text(.15,.81,['Miss distance = ' num2str(Rtm(mi)) ' ft'],'sc’); 
text(. 15,-77, ['Quant Error = ' nuni2str(QE) ' ft'],'sc'); 
text(.15,.73,['Inter Error = ' nuni2str(r) ' ft'], 'sc'); 
text(.15,.69,['Sigma(It) = ' num2str(ub(It)* 180/pi) ' deg’],'sc'); 
text(.15,.65,['Max Acc. Limit = ' num2str(maxac) ' ft/sec^2'],'sc'); 
title('Control: Proportional Navigation') 
xlabel('Direction 1 (ft)'),ylabel('Direction 2 (ft)'); 
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5. BB3RD0RD.M 



% BB3RDORD.M 3/7/91 

% This is a simulation of a 3rd order system, forward time. 

% This version allows varied N values. 

% written by Colin R. Cooper 



xl0=-.5; % Define initial condirions 

x20=0; 

x30=0; 



N=l; % Set the magnitude of the control effort. 

Tf=2.5; % Set maximum time of simulation. 

dt=.002; % Set the simulation step size. 



A=[0 1 0;0 0 1;0 0 0]; % Define the State Equations 

B=[0 0 1]'; 

C=[l 0 0]; 



[Phi,Del]=c2d(A,B,dt); 
kmax=Tf/dt +1; 



% Discretize the system. 

% Max integer value for the simulation. 



% Prepare storage vectors. 

u=zeros(l,kmax); 

x=zeros(3,kmax); 

y=zeros(l,kmax); 

w=zeros(l,kmax); 

f=zeros(l,kmax); 

time=zeros( 1 ,kmax) ; 
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x(:,l)=[xl0;x20;x30]; 



% Define initial conditions in state vectors 



% Begin loop for simulation, 
for (i=l:kmax-l) 

w(i)=sign(x(2,i)+x(3,i)*abs(x(3,i))/(2*N)); % Defining the switching law. 

f(i)=x(2,i)+w(i)*(x(3,i)'^2)/(2*N); 

u(l,i)=-N*sign(x(l,i) + (x(3,i)'^3)/(3*N'^2) + w(i)*x(3,i)*x(2,i)/N +... 
f(i)*abs(f(i)/N)^.5); 

x(:,i+l) = Phi*x(:,i) +Del*u(l,i); % Calculate the state values. 

y(l,i+l) = C*x(:,i+l); 

time(i+l)= time(i) + dt; % Store time vector, 

end; 

% Plot the switching law and its components, 
clg, axis([0 2.5 -1.5 1.5]); 

plot(time,u, time, .75*w, time, f) % w is scaled to distinquish it from u. 

title('Control Laws vs Time') 

pause 

% Plot the 3-Dimensional view of the simulation from 45 deg. azimuth 
% and 45 deg. elevation angle. (Pos. xl vector is out of the screen 
% and towards the left, Pos. x2 is out of the screen and towards right, 
clg, axis; 

plot3d(x(l,:),x(2,:),x(3, 0,45,45); 
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6. ADJOINT.M 



% ADJOINT.M is a simulation of a controlled system using formard time simulation 
% and reverse time or adjoint solution. 



% written by Colin R. Cooper 

A = [0 1 0;0 0 1;0 0 0]; 

B = [0 0 1]’; 

C = [l 0 0]; 
xlO = 2; 
x20 = 0; 

N = 1 ; % Maximum control effort. 

Tf = 4; 
dt = .005; 

tsf = x20/N + sqrt(xl0/N + .5*(x20/Nr2); 
tsa = sqrt(xl0/N + .5*(x20/N)^2); 
kmax = Tf/dt + 1; 

% Define the Storage Vectors 
X = zeros(3,kmax); 
y = zeros(l,kmax); 
time = zeros(l,kmax); 
imp = zeros(l,kmax); 
imp(l) = -N/dt; 

imp(round(tsf/dt)+l) =2*N/dt; 
x(l:2,l) = [xl0;x20]; 
[Phi,del]=c2d(A,B,dt); 

for i = 1 :kmax - 1 

x(:,i+l) = Phi*x(:,i) + del*imp(i); 
y(l,i+l) = C*x(:,i+l); 



11 Mar. 1991 

% Define Forward Time State Equations 

% xl initial condition. 

% x2 initial condition. 

% Length of simulation. 

% Time increment for simulation. 

% Forward Time Switching Time 
% Adjoint System Switching Time 
% Max length of storage vectors. 

% States of the system 
% Output state vector 

% Impulse vector for control times. 

% First pulse at t = 0 
% Second pulse at t = tsf 
% Set initial conditions for x. 

% Discretize the state equations 

% Begin the simulation loop 



72 



time(i+l) = time(i)+<it; 



end 

miss = min(sqrt(x(l,:).^2 + x(2,:).'^2)); % Find the min. miss distance 

tff = time(find(sqrt(x(l,:).'^2 + x(2,:).'^2)== ... % Find the dme of the min. miss 

min(sqrt(x(l,:).''2 + x(2,:).^2)))); 

clg,plot(x(l,:),x(2,:),'-w’),grid % Plot the Phase Plane 

title('Forward Time Phase plot') % Lable the graph and display the desired 

xlabel(’Xr),ylabel('X2') % information 

text(.6,.80,[’min miss = ' num2str(miss)],'sc’) 
text(.6,.77,['tsf = ' num2str(ts0 ’ sec.’], 'sc') 
text(.6,.74,['tff = ’ num2str(tff) ' sec.'],'sc') 
pause 

% Now for the Adjoint System, 
xa = zeros(3,kmax); 
time = zeros(l,kmax); 
impa = zeros(l,kmax); 
impa(l) = N/dt; 

impa(round(tsa/dt)+l) = -2* N/dt; 

[Phi,del]=c2d(A’,C’,dt); 

for i = 1 :kmax - 1 

xa(:,i+l) = Phi*xa(:,i) + del*impa(i); 
ya(l,i+l) = B’*xa(:,i+1); 
time(i+l) = time(i) + dt; 
end; 

missa = min(sqrt((xa(3,:)-xl0).'^2 + (xa(2,:)-x20).'’^2)); % Find the min. miss distance 

% from the initial conditions 

tfa = time(find(sqrt((xa(3,:)-xl0).'^2 + (xa(2,:)-x20).'"'2)== ... % Find the time of the 

min(sqrt((xa(3,:)-xl0).^2 + (xa(2,:)-x20).'''2)))); % min miss distance 



% Define storage vectors 

% First impulse 
% Second impulse 
% Discretize the Adjoint System 

% Begin the simulation loop 



plot(xa(3,:),-xa(2,:),'-w'),grid 
title('Adjoint Solution Phase plot') 
xlabel(’X3'),ylabel('X2’) 
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% Plot the Phase Plane 
% Label the graph and display the desired 
% information. 



text(.6,.8,['min miss = ' num2str(missa)],'sc') 
text(.6,.77,['tsa = ' num2str(tsa) ' sec.'], 'sc') 
text(.6,.74,['tfa = ' num2str(tfa) ' sec.'],'sc') 
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7. ADJSIM.M 



% ADJSIM.M Adjoint solution to missile intercept problem. 



% written by Colin R. Cooper 


11 Mar. 1991 


n = 4; 


% Proportionality Constant 


V = 5000; 


% Closing Velocity 


B = [0 000 0 1]'; 


% Define State Matrices 


C = [-l 000 1 0]; 




Tf= 1; 


% Set first value for incremented times 


dt = .01; 




% This outer loop runs the complete forward time system for each value of Tf, storing 


% the output into vectors at the end of the loop. 


for J = 1:4 




kmax = Tf/dt + 1; 


% Maximum index for storage vectors 


imp = zeros(l,kmax); 


% Create vector for impulsive input 


imp(l) = 1/dt; 


% Define the pulse at t = 0 


X = zeros(6,kmax); 


% Create storage vectors for states 


y = zeros(l,kmax); 




time = zeros(l,kmax); 




% Forward Time Simulation. 




count = 0; 


% Counter to indicate the computer is busy 


for i = 1 :kmax - 1 


% Begin simulation loop 


count = count + 1; 




ki = l/(v*(Tf - time(i) + le-12)); 


% 1/Range 


A = [0 1 000 0 


% Define the time varying A matrix 


0 0n*v00 0 




000100 
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-100*ki 0-100 -20 100*ki0 
000001 
0 0 00 00 ]; 

Phi = eye(6) -i- A*dt ,5*A^2*dt^2; 

x(:,i-t-l) = Phi*x(:,i) + B*dt*imp(i); 
y(:,i+l) = C*x(:,i+l); 
time(i+l) = time(i) + dt; 
if count == 50 

disp('working'); 
count = 0; 
end 
end 

Tend(j)=Tf; 
yendO) = y(length(y)); 
eval(['y’,num2str0),’ = y;'J); 
eval(['t',num2str(j),' = time;’]);l 
Tf = Tf-i- 1 
end 

% Adjoint Simulation. 

Tf=4.25; 
kmax = Tf/dt -i- 1; 
imp = zeros(l,kmax); 
imp(l) = 1/dt; 
xa = zeros(6,kmax); 
ya = zeros(l,kmax); 
time = zeros(l,kmax); 
count = 0; 
for i = l:kmax - 1 
count = count -i- 1; 
ki = l/(v*(time(i) + le-12)); 

A = [0 1 0 0 0 0 



% Discretize the matrix with 2nd order 
% expansion - drop higher terms 



% Counter indicates computer is busy 



% End simulation loop 
% Record final time 

% Record corresponding final value ymiss 
% Save output vector 
% Save corresponding time vector 
% Increment to next Tf 
% End outer forward time loop 

% Set a simulation time 
% Maximum index for vectors 
% Create vector for impulsive input 
% Define the impulse at t = 0 
% Create storage vectors 



% Begin simulation loop 
% 1/Range 

% Define the time varying A matrix 



76 



00 n*vOOO 
000100 

-100*ki 0-100 -20 100*ki0 
000001 
00 00 00 ]; 

Phi = eye(6) + A’*dt + .5*A’'^2*dt'^2; % Discretize the matrix 

xa(:,i+l) = Phi*xa(:,i) + C*dt*imp(i); 
ya(:,i+l) = B'*xa(:,i+D; 
time(i+l) = time(i) + dt; 
if count == 50 
disp('working'); 
count = 0; 
end 

end % End adjoint simulation loop 



plot(tl,yl,’-w',t2,y2,'-w',t3,y3,'-w',t4,y4,'-w',time,ya,'-w') % plot output 

xlabel(Time (sec)'),ylabel('Ymiss (ft)'), grid 
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8. INTERP.M 



function r = interp(xm,xt) 

% INTERP.M will return an interpolated value for the min.miss distance given the states 
% for the two intercept values. 

% Written by Colin R. Cooper 26 Mar 1991 

% Increase the sample rate by a factor of 100 in crossover region. 

dax = (xm(3,2)-xm(3,l))/100; 

day = (xm(l,2)-xm(l,l))/100; 

dbx = (xt(3,2)-xt(3,l))/100; 

dby = (xt(l,2)-xt(l,l))/100; 

% Define the storage vectors for 100 data points, 
ax = zeros(l,100); 
ay = zeros(l,100); 
bx = zeros(l,100); 
by = zeros(l,100); 

% Set initial values for each vector. 

ax(l) = xm(3,l); 

ay(l) = xm(l,l); 

bx(l) = xt(3,l); 

by(l) = xt(U); 

% Assuming a straight line trajectory with constant velocity in the crossover region, create 
% the interpolation data sets, 
for i = 1:99 

ax(i+l) = ax(i) + dax; 
ay(i+l) = ay(i) +day; 
bx(i+l) = bx(i) + dbx; 
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by(i+l) = by(i) + dby; 



end 



% Find the closest point of approach of the two line segments, 
r = min(sqrt((ax - bx).'^2 + (ay - by).^2)); 
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9. PL0T3D.M 



% PLOT3D is a 3d plotting function allowing rotation and elevation adjustments. The plot 
% shows a 3-D curve and its projection onto the X-Y Plane. 

% X,Y,and Z data must be passed, and if no azimuth or elevation values are passed 
% they default to AZ = 45°, EL = 30°. The Azimuth is the angle of rotation of the view 
% angle about the Z-Axis. AZ = 0° is looking straight down the X-Axis at the Y-Z Plane. 
% The elevation is the angle from which the plot is viewed. EL = 90° is looking down the 
% Z-Axis at the X-Y Plane. The elevation angle can vary from -90° to 90°. 

% The tick marks on the axis will default to 10 marks per axis and the values will be 
% displayed on the screen. To define the values of the tick marks a three element vector 
% must be passed containing [dx dy dz]; 

% Axis values will be calculated and fixed unless an axis vector is passed to the 
% program: [xmin xmax ymin ymax]. 

% The Transformed 2-D vectors V and Vs are returned, where V is the 3-D curve and 
% Vs is the Projection on the X-Y Plane. 

% 

% Example : [V,Vs] = plot3d(x,y,z,-45,30,dx,Ax) 

% written by Colin Cooper 4/26/91 

function [V,Vs]=plot3d(x,y,z,az,el,dx,Ax) 

if nargin < 5, el = 30; end 
if nargin < 4, az = 45; end 



az=-az*pi/180; el=el*pi/180; 

alpha = 45*pi/180; beta=30*pi/180; % Angles for mapping onto 2D 

% alpha = rot. Beta = elv. 

Sl=[sin(az) cos(az) 0 

-sin(el)*cos(az) sin(el)*sin(az) cos(el)]; 

S2=[sin(az) cos(az) 0 

-sin(el)*cos(az) sin(el)*sin(az) 0]; % Trans for Projection. 
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[r,c]=size(x); 
ifr> 1 

X = x'; y = y'; z = z'; 
end 

V = Sl*[x;y;z]; 

Vs = S2*[x;y;z]; 

if nargin < 6 

dx(l) = (max(x)-min(x))/10; 
dx(2) = (max(y)-min(y))/10; 
dx(3) = (max(z)-min(z))/10; 
end 

axl = [min(x) max(x); 0 0; 0 0]; 
ax2 = [0 0; min(y) max(y); 0 0]; 
ax3 = [0 0; 0 0; min(z) max(z)]; 

axlt = [fliplr(0;-dx(l):niin(x)) 0;dx(l):max(x) 
zeros([0:-dx(l):min(x) 0:dx(l):max(x)]) 
zeros([0:-dx(l):min(x) 0:dx(l):max(x)])]; 
ax2t = [zeros([0:-dx(2):min(y) 0:dx(2):max(y)]) 
fliplr(0:-dx(2);min(y)) 0:dx(2);max(y) 
zeros([0:-dx(2):min(y) 0:dx(2):max(y)])]; 
ax3t = [zeros([0;-dx(3):min(z) 0:dx(3):max(z)]) 
zeros([0:-dx(3):min(z) 0:dx(3):max(z)]) 
fliplr(0:-dx(3):min(z)) 0:dx(3):max(z)]; 

axl = Sl*axl; 
ax2 = Sl*ax2; 
ax3 = Sl*ax3; 
axlt = Sl*axlt; 
ax2t = Sl*ax2t; 
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ax3t = Sl*ax3t; 



minx = min([V(l,:) Vs(l,:) axl(l,:) ax2(l,:) ax3(l,:)]); 
if minx == 0, minx = -1; end 

maxx = max([V(l,:) Vs(l,:) axl(l,;) ax2(l,;) ax3(l,:)]); 
if maxx == 0, maxx = 1; end 

miny = min([V(2,:) Vs(2,:) axl(2,;) ax2(2,:) ax3(2,:)]); 
if miny == 0, miny = - 1 ; end 

maxy = max([V(2,;) Vs(2,:) ax 1(2,:) ax2(2,:) ax3(2,:)]); 
if maxy == 0, maxy = 1; end 



cig 

axis('square'); 
if nargin < 7, 

axis([minx-.3*abs(minx) maxx+.3*abs(maxx) miny-.3*abs(miny) . 
maxy+.3*abs(maxy)]); 

else 

axis(Ax); 

end 

b = axis; 
hold on 

plot(axl(l,:),axl(2,:),'-w',axlt(l,:),axlt(2,:),'+w',... 

ax2(l,:),ax2(2,:),'-w',ax2t(l,:),ax2t(2,:),’+w',... 

ax3(l,:),ax3(2,;),'-w',ax3t(l,:),ax3t(2,:);+w') 

Iv = length(V); 

plot(V(l,:),V(2,:),'-w',Vs(l,:),Vs(2,:),'-w',... 

[V(l,l:15:lv);Vs(l,l:15:lv)],[V(2,l:15:lv);Vs(2,l:15:lv)],':w') 
plot([b(l) b(2) b(2) b(l) b(l)],[b(3) b(3) b(4) b(4) b(3)],'-w'); 
hold off 

text(.22,.08,['AZ = ’ num2str(-az* 180/pi) ’°’],’sc'); 
text(.72,.08,['EL = ' num2str(el* 180/pi) '“'l/sc'); 
text(.22,.88,['dx = ' num2str(dx(l)) ],'sc'); 
text(.50,.88,['dy = ' num2str(dx(2)) ],'sc'); 
text(.72,.88,['dz = ’ num2str(dx(3)) ],'sc'); 
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